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Abstract. Recent cosmological data favour additional relativistic degrees of freedom beyond 
the three active neutrinos and photons, often referred to as "dark" radiation. Light sterile 
neutrinos is one of the prime candidates for such additional radiation. However, constraints 
on sterile neutrinos based on the current cosmological data have been derived using simplified 
assumptions about thermalisation of v s at the Big Bang Nucleosynthesis (BBN) epoch. These 
assumptions are not necessarily justified and here we solve the full quantum kinetic equations 
in the (1 active + 1 sterile) scenario and derive the number of thermalised species just before 
BBN begins (T ~ 1 MeV) for null (L = 0) and large (L = 10~ 2 ) initial lepton asymmetry 
and for a range of possible mass-mixing parameters. We find that the full thermalisation 
assumption during the BBN epoch is justified for initial small lepton asymmetry only. Partial 
or null thermalisation occurs when the initial lepton asymmetry is large. 
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1 Introduction 

Sterile neutrinos are hypothetical SU(2) x U(l) singlets. They are supposed to mix with one 
or more of the active states without interacting with any other particle. Low-mass sterile 
neutrinos have been invoked to explain the excess v e events in the LSND experiment [1-3] as 
well as the MiniBooNE excess events in both neutrino and antineutrino channels. Interpreted 
in terms of flavour oscillations, the MiniBooNE data require CP violation and thus no less 
than two sterile families [4-6] or additional ingredients such as non-standard interactions [7] . 
Recently a new analysis of reactor v e spectra and their distance and energy variation [8-10] 
suggested indication for the possible existence of eV-mass sterile neutrinos. However the 
IceCube collaboration excluded part of the parameter space [11]. 

The most recent analysis of cosmological data suggest a trend towards the existence of 
"dark radiation," radiation in excess with respect to the three neutrino families and pho- 
tons [12-14]. The cosmic radiation content is usually expressed in terms of the effective 
number of thermally excited neutrino species, N e g. Its standard value, N e g = 3.046, slightly 
exceeds 3 because of e + e~ annihilation providing residual neutrino heating [15]. The Wilkin- 
son Microwave Anisotropy Probe (WMAP) collaboration found JV e g- = 4.34^q|| based on 
their 7-year data release and additional LSS data [16] at la. Including the Sloan Digital Sky 
Survey (SDSS) data release 7 (DR7) halo power spectrum, [13] found N e g = 4.78^-^;^ at 
2a. Measurements of the CMB anisotropy on smaller scales by the ACT [17] and SPT [18] 
collaborations also find tentative evidence for a value of Neg higher than predicted by the 
standard model (see also [19-23] for recent discussions of N e g). 

Also, cosmological constraints coming from big bang nucleosynthesis (BBN) suggest 
that the relatively high He abundance can be interpreted in terms of additional radiation 
during the BBN epoch [24, 25]. Low-mass sterile neutrinos have been considered among pos- 
sible candidates for the extra-radiation content [26-28]. The cosmic microwave background 
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anisotropies and big-bang nucleosynthesis in combination seem to favor an excess of radia- 
tion compatible with one family of sub-eV sterile neutrinos [14, 27-29]. On the other hand, 
eV-mass sterile neutrinos are cosmologically viable only if additional ingredients are included 
since otherwise sterile neutrinos would contribute too much hot dark matter [26] (see also 



However, cosmological constraints during the BBN epoch have usually been derived 
under the assumption that the extra sterile neutrino families were fully thermalised [26]. 
However, the validity of this assumption is not a priori clear and some preliminary studies 
[31, 32] already pointed toward this direction. It was shown in [33] that for plausible values 
of the mass and mixing parameters, and initial lepton asymmetries not excluded by current 
observations there are cases where little or no thermalisation occurs. For the charged fermions 
of the standard model the particle anti-particle asymmetry is known to be of order 10 -10 . 
For neutrinos, however, no such bound exists, and the asymmetry can be many orders of 
magnitude larger without violating observational constraints. In the standard model with no 
sterile states the upper bound on the neutrino chemical potential is of order fi/T < few x 10 -2 
[34-39] and while no exact bound has been derived in models with sterile neutrinos, we expect 
that the upper bound is of the same order of magnitude. 

The purpose of this paper is to quantitatively derive the amount of thermalisation as a 
function of neutrino parameters (mass, mixing, and initial lepton asymmetry). We solve the 
full quantum kinetic equations in the 1 active+1 sterile approximation, calculate the effective 
number of thermalised species just before BBN starts (at T ~ 1 MeV) and define under which 
conditions the thermalisation hypothesis holds. The assumption of (1 + 1) families to evaluate 
the thermalisation degree is justified for small lepton asymmetries since the resonances in the 
active sector are decoupled from the conversions occurring in the active-sterile sector due to 
the larger mass difference. However, for large asymmetries active-sterile conversion is delayed 
and can occur simultaneously with active-active conversion. While this does not qualitatively 
change the overall picture there are some issues which we will return to in Section 3. 

In our study, we calculate the number of thermalised extra families for the allowed mass- 
mixing parameter space for different initial lepton asymmetries. In Section 2 we introduce 
the adopted formalism and the quantum kinetic equations. In Section 3, we present our 
results for initial null and large (L = 10~ 2 ) lepton asymmetry. Conclusions and perspectives 
are presented in Section 4. 

2 Equations of motion 

In this section, we introduce the quantum kinetic equations (QKEs) governing the evolution 
of neutrinos in the early universe [40-45] . We adopt a mapping of the Bloch vectors in terms 
of new vectors related to the active and sterile species grouping large and small dynamical 
variables. 

2.1 Quantum Kinetic Equations 

We consider oscillations of one active flavour v a (with a = e or fi, r) with a sterile neutrino 
state v s . Denoting with 6 S the mixing angle in vacuum and with u\ and V2 the two mass 
eigenstates, separated by the mass difference 5m 2 s , we have: 



[30]). 



u a = cos9 s v\ — sin S v<i , 
v s = sin^ s z^i + cos Q S V2 ■ 




(2.2) 
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In what follows we will refer 5mj > as the normal hierarchy scenario (NH) and 5mj < 
as the inverted hierarchy scenario (IH). Structure formation data strongly disfavour models 
with a total thermalised neutrino mass (the sum of all fully thermalised mass states) in 
excess of 0.5-1 eV. Given that all the active states are fully thermalised this disfavours the 
inverted hierarchy for sterile masses above 0.2-0.3 eV. However, for masses below this the 
inverted hierarchy is not disfavoured and for completeness we study the same mass and 
mixing parameter space for both NH and IH. 

In order to describe the evolution of sterile neutrinos in the early universe, we use the 
density matrix formalism and we express the density matrix associated with each momentum 
p in terms of the Bloch vector components (Po,P) = (Po, P x , P y , P z ) [40, 41, 43], 

p=^fo(P + P-v) , p = \fo{P~o + P • a) , (2.3) 

where a are the Pauli matrices and /o = l/(l + e P//T ) is the Fermi-Dirac distribution function 
with no chemical potential. The neutrino kinetic equations in terms of the components of 
the Bloch vectors for each momentum mode are: 

P = V x P - D(P x x + P y y) + P z , (2.4) 



Po = r 



f -f~\{Po + P z 
JO 2 



(2.5) 



where the dot denotes the time derivative (dj = dt — Hpd p , with H the Hubble parameter) 
and / eq = 1/(1 + e^-^l T ). 

Defining the comoving momentum x = p/T, the vector V has the following components 

V x = ^sm29 s , (2.6) 

V y = 0, (2.7) 
V z = V + V 1 + V L . (2.8) 



and 



Am 

Vo = -2^T COS29s > (2 - 9) 

v ^ = -^m xT5[n - +n ^ (2 - 10) 



v L 



^^±G F T 3 L( a \ (2.11) 



Here, g^ T = 1 for v^ T -v s mixing, g e = 1 + 4 sec 2 9w/( n u e + nv e ) for u e -v s mixing and 9w is 
the Weinberg angle. The dimensionless number densities n v rp a \ are the equilibrium active 
neutrino (antineutrino) densities normalised to unity in thermal equilibrium. The effective 
neutrino asymmetries are defined by 

L [e) = Q + 2 sin 2 9w\ L e + Q - 2 sin 2 9 W \ L p - h n + 2L Ue + + L Ut , (2.12) 

Z>) = L( £ ) -L e - L Ve +L, M , (2.13) 
£ W = L {e) - Le - L Ve + L Vt , (2.14) 
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where Lf = (rif — nj)Nf/N^ with Nf (A r 7 ) the integrated active (photon) number density 
in thermal equilibrium. The potential Vl, defined as in Eq. (2.11), is the leading order 
contribution to V z . The V\ term is the finite temperature correction and for example in the 
case of v e —v s mixing it includes coherent interactions of v e with the medium through which 
it propagates. The condition for a matter induced resonance to occur is V z = 0, and because 
V z depends on any non-zero lepton asymmetry can have dramatic consequences for 
oscillation driven active-sterile neutrino conversion. In Appendix A we discuss the location 
of resonances in detail for all possible values of mass, mixing, and lepton asymmetry. 

A detailed derivation of the quantum kinetic equations is presented in [42, 46]. Here 
we choose to adopt minimal assumptions on the collision terms. In particular, the term D is 
the damping term, quantifying the loss of quantum coherence due to v a collisions with the 
background medium. For example, considering u e , the elastic contribution should come from 
the elastic scattering of v e with e~ and e + and with the other active flavours v a and v a . The 
inelastic contribution comes from the scattering of u e with v e (producing e~ and e + or v a and 
Va). In terms of the Bloch vectors such terms have the effect of suppressing the off-diagonal 
elements of the density matrix (P x ,y)- The effective potentials contributing to this term have 
been previously calculated [46-48] and if thermal equilibrium is aasumed and the electron 
mass neglected, it is approximately half the corresponding scattering rate T [40, 42, 49] 



The evolution of Pq is determined by processes that deplete or enhance the abundance 
of v a with the same momentum and its rate of change receives no contribution from coherent 
v a -v s oscillations. The repopulation term r(/ eq //o — 1/2(Pq + Pz)) is an approximation for 
the correct elastic collision integral [49] with 



where C e ~ 1.27 and Cut — 0.92 [41]. Note that the term including the effective collision rate, 
r, is an approximation to the full momentum dependent scattering kernel which repopulates 
neutrinos from the background plasma. The full expression has been derived in [42] . In [49] it 
was proven that the general form of D (and T) exactly reduces to Eqs. (2.15,2.16) for weakly 
interacting species in thermal equilibrium with zero chemical potential, and that it is the zero 
order approximation for particles with non-null chemical potential. The respective equations 
of motion for anti-neutrinos can be found by substituting = —L^ and fi = —\i in the 
above equations. In our treatment we have not included the rate equations for the electrons 
and positrons since we are assuming that all the species electromagnetically interacting are 
kept in equilibrium. 

2.2 Mapping with the active and the sterile variables 

We can distinguish among large and small linear combinations of the dynamical variables in 
the particle and antiparticle sector to simplify the numerical treatment. For each momentum 
mode, we define for each component i (with i = 0, x, y, z) of the Bloch vector 




(2.15) 



r = C a GpxT 5 





Pi±Pi- 



(2.17) 
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We also separate active (a) and sterile (s) sectors 



P: 



P d 



fo 
pt 
fo 



Therefore, in terms of the new vectors Eqs. (2.4, 2.5) become 
Pa =V x P± + T[2ff q /f Q -Pt] , 



s 



-V P ± 

V X fy 



(V + Vi)P±-V L P?-DP, 



± 



(V + V^Pt + V L PT - - Vx {Pt - Pt) ~ DP, 



(2.18) 
(2.19) 

(2.20) 
(2.21) 
(2.22) 

(2.23) 



P 

where we have defined /± = f cq (p, n) ± f cq (p, -/i). 

The lepton number can be directly calculated from the integral over the difference 
between the neutrino and the antineutrino distribution functions, i.e. P~: 



L («) 



8C(3) 



dxx p~ 



1 



8C(3) 



dxx /o^a • 



(2.24) 



However, since the repopulation term is approximated by Eq. (2.16) which does not 
explicitly conserve lepton number we independently evolve as in [50] using an evolution 
equation where the repopulation term does not enter. Taking the time derivative of Eq. (2.24) 
and ignoring the repopulation part of P~ , the evolution equation for is 



£(«) 



8C(3) 



dxx 2 f V x p- 



(2.25) 



Note that, in kinetic equilibrium, fj,, or rather the degeneracy parameter £ = p/T, is 
related to the lepton number through the integral over [34] 



L (a) 



4C(3) 



dx x 



1 



1 



1 + e x -Z 1 + e x +Z 



12C(3) 



(7r 2 £ + e 3 ). (2.26) 



This is a third order equation, and using Chebyshev's cubic root, one can extract the corre- 
sponding and expression for £ valid for any using trigonometric functions: 



-2vr 

7T 



sinh [ — arcsinh 



W3C(3) L(a) 



7T" 



(2.27) 



In order to numerically solve the QKEs, we define the momentum grid in comoving 
coordinates (x = p/T). Therefore the grid becomes stationary and the partial differential 
equations become ordinary differential equations coupled through integrated quantities only. 
Using the temperature T as the evolution parameter, time derivatives, dt, are replaced by 
— > —HTOt in the above equations, provided that the time derivative of the effective number 
of degrees of freedom can be ignored. 
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3 Results: thermalised sterile species 



The fraction of sterile thermalised species is defined as 

fdxx 3 f P+ , . 

= ijW • (31) 

However, the total amount of radiation is given by the sum of active and sterile energy 
densities 

f dxx 3 f (P+ +P+-4) . . 

5N " = 1 — wif • (3 - 2) 

Note that when the active state is in thermal equilibrium (P+ = 4), <5iV e ff s = 5N c q. When 
is large, the sterile sector may be populated so late that the active sector does not have 
time to repopulate before it decouples. In this section, we discuss the fraction of thermalised 
species for initial = and = 10 -2 and for a range of (5m 2 , sin 2 26 s ). 

In terms of late-time cosmological constraints on light neutrinos both SN e Q s and 5N c g 
can be relevant quantities. Models with a modified light neutrino sector are most often 
parametrised in terms of the neutrino mass, m„, and N c g in such a way that N e g neutrino 
species all share the same common mass m v (i.e. it is assumed that the mass spectrum is 
degenerate). However, in models with a single sterile state one instead has either 8N e g B 
steriles with mass m s and 3.046 + SN e g — 5N c g s massless active states (NH) or 3.046 + 
8N e ff — 5N e fi jS massive active state with degenerate mass and 5N e s tS massless sterile states 
(IH). These two cases are different when it comes to structure formation and should in 
principle be treated separately (see e.g. [51] for a discussion about this point). Since the goal 
of this paper is to calculate 5N e g, not to provide quantitative constraints on specific models, 
we simply use 5N e g from this point on. 

3.1 Numerical solution of the quantum kinetic equations 

Solving the quantum kinetic equations numerically is non-trivial task. The number of differ- 
ential equations are roughly 8iV where N is the number of momentum bins, and since the 
resonances can be very narrow we need a few hundred points to obtain good precision. There 
are many vastly separated time-scales involved, so the problem is stiff, and once L changes, 
the system becomes extremely non-linear. We used two different solvers, one based on the 
numerical differentiation formulae of order 1 — 5 (ndf 15) due to Shampine [52], and one based 
on the fifth order implicit Runge-Kutta method RADAU5 due to Hairer and Wanner [53]. If 
the maximum order of the first method is reduced to two, both solvers are L-stable, and thus 
excellent for stiff problems. Because of the large number of equations and the sparsity of the 
Jacobian we must use sparse matrix methods for the linear algebra operations needed in both 
solvers. For this purpose, we are employing a small sparse matrix package based on [54]. 

To sample the momentum-space in an optimal way we are mapping the x-interval 
[xmm, ^rnax] to a u- interval [0; 1] by 

/ \ X m ax -^ext / \ 

U[X) = X , (3.3) 

3-max ^min % ^ext 

where x ex t is the extremal point of some moment of the Fermi-Dirac distribution. We chose 
the values x m - m = 10 -4 , x ex t = 3.1 and x max = 100, and then sampled u uniformly. This is 
the same mapping employed in [55], but they go one step further and introduce an adaptive 
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grid that follows the resonances. This is not necessary for this project since our mixing angles 
are comparably larger, and we are not looking at chaotic amplification of an initially small 
value of L. 

We evolved the system from an initial temperature of 60 MeV to a final temperature of 
1 MeV for the following grid of masses and mixing angles: 

lfT 3 eV 2 < 5m 2 s < 10 eV 2 and 10~ 4 < sin 2 26 s < 10" 1 for = , (3.4a) 

10" 1 eV 2 < 5m 2 s < 10 eV 2 and 10~ 3 ' 3 < sin 2 26 s < 10" 1 for L (a) = 10~ 2 . (3.4b) 

We ran the complete grids for different number of momentum bins, different accuracy pa- 
rameters and both differential equation solvers with no noticeable difference. 

3.2 Sterile neutrino production for zero lepton asymmetry 

The simplest case, and the one most often studied in the literature, is the one where the 
lepton asymmetry is zero. For LW = 0, the evolution of P- + is decoupled from [sec 
Eqs. (2.20,2.23)] and the asymmetry remains zero for the whole evolution (as can be seen 
from Eq. (2.25)). 

From Eqs. (A. 4, A. 8) in Appendix A it can be seen that there is either no resonance 
(NH) or that the resonances are identical for neutrinos and anti- neutrinos (IH). As it is well 
known, in IH the resonance propagates to higher values of x as the universe expands and 
eventually covers the entire momentum distribution of neutrinos. 

In Fig. 1 we show the fraction of thermalised neutrinos, 5N c h, for the range of mixing 
parameters given in Eq. (3.4a) with initial asymmetry = 0. The top panel shows the 
normal hierarchy, 5m? s > 0, and the bottom panel the inverted hierarchy, <5m 2 < 0. The 
smaller parameter space described by (3.4b) is denoted with a dashed rectangle to facilitate 
comparison with the results presented in Sec 3.3. 

We mark with a green hexagon the best fit point of the 3 + 1 global analysis presented 
in [56], obtained from a joint analysis of Solar, reactor, and short-baseline neutrino oscillation 
data (<5m 2 , sin 2 29 s ) = (0.9 eV 2 , 0.089). For that point 5N D g = 1 in both hierarchies, i.e. 
complete thermalization occurs. In addition we show the parameter range preferred by CMB 
and large scale structure (LSS) data. The 1 — 2 — 3a contours have been obtained interpolating 
the likelihood function obtained in [27] for each fixed 5m 2 and N e g. In both cases the lower 
left corners of parameter space where little thermalization occurs are disfavoured because of 
the CMB+LSS preference for extra energy density. 

It is also of interest to see how the thermalization proceeds as a function of temperature. 
In Fig. 2 we show the evolution of 5N e g as a function of temperature for the NH scenario for 
a variety of different <5m 2 and sin 2 2# s . For the non-resonant NH, the thermalization rate of 
sterile neutrinos is approximately T s ~ | sin 2 29 S T. The maximum thermalisation rate occurs 
at a temperature of approximately 10 (5m 2 ) 1 / 6 MeV and the final 8N e s depends only 

on sin 2 26 m at that temperature (see [41] for a detailed discussion). In the top panel of Fig. 2 
this behaviour can be seen. For very large vacuum mixing T s /H > 1 already before T max 
such that complete thermalisation has occurred already before T max reached. For smaller 
mixing T s /H never exceeds 1 and even though thermalisation proceeds fastest around T max 
it is never fast enough to equilibrate the sterile states. 

In the bottom panel the change in T max as <5m 2 varies is evident, and provided that 
Tmax is higher than the active neutrino decoupling temperature the vacuum mixing in this 
case is large enough that complete thermalisation always occurs. For the non-resonant case 
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Figure 1. lso-5N c g contours in the sin 20 s 



plane for LW =0 and 8m 2 s > (top panel) and 



5m a < (bottom panel). The green hexagon denotes the v s best-fit mixing parameters as in the 3 + 1 
global fit in [56]: (6m 2 s , sin 2 20 s ) = (0.9 cV 2 , 0.089). The 1 - 2 - 3cr contours denote the CMB+LSS 
allowed regions for v s with sub-eV mass as in [27]. In order to facilitate the comparison with the 
results presented in Sec 3.3, a dashed rectangle denotes the parameter-space described by (3.4b). 



the end result is that isocontours of 5N e g always lie at constant values of 5m 2 s sin 4 20 s , as can 
be seen in the top panel of Fig. 1. 

In the inverted hierarchy the resonance conditions are always satisfied. Therefore, we 
expect full thermalization for a larger region of the mass-mixing parameters than in NH, as 
confirmed in Fig. 1. In this case, thermalisation may proceed through resonant conversions 
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T (MeV) 

Figure 2. Top panel: 5N c g as a function of the temperature for four different mixing angles 
(sin 2 26 s = 1(T 4 ,2 x 1CT 3 ,5 x lO^lCT 1 ) and fixed mass difference (<5m 2 = 0.93 eV 2 ). Bot- 
tom panel: SN c g as a function of the temperature for four different mass differences (5m 2 . = 
1CT 3 , 3.5 x 10~ 2 , 9.3 x lfr\ 10 eV 2 ) and fixed mixing angle (sin 2 26 s = 0.051). Thermalisation begins 
earlier and is more effective for larger mass differences and for larger mixing angles. 

alone. For illustration, we choose the point of Fig. 1 with (5m 2 s , sin 2 6 S ) = (-3.3 eV 2 , 6x 10~ 4 ) 
for which 8N e s = 0.55 and we show the percentage of active (-/V a ) and sterile (N s ) neutrinos 
as a function of x for different T in Fig. 3. The thermalisation is not complete and it is nearly 
instantaneous as the resonance moves through the momentum spectrum and the resulting 
dip in the active sector is quickly repopulated from the background. 

We have presented results for = only, but the case of = shows exactly the 
same trend as in Fig. 1. However, the region with <5iV e ff = 1 is slightly smaller than the one 
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Figure 3. Temperature evolution of active and sterile neutrino distributions for the resonant case 
(<5m 2 , sin 2 9 a ) = (-3.3 eV 2 , 6 x 10" 4 ) and L<» = 0. 



shown in Fig. 1. This is due to the fact that v^s have a larger potential than (because of 
the charged current interaction contribution) and therefore resonances occur at slight lower 
temperatures. 

3.3 The case of large initial lepton asymmetry 

We now discuss the thermalisation degree for initial large lepton asymmetry. In principle, 
one would expect a lepton asymmetry of the same order of magnitude as the baryon asym- 
metry (77 ~ 1CP 10 ). However, since neutrinos are neutral particles, = 10~ 2 - 10 _1 is 
not presently excluded [36, 37, 57] by the requirement of charge neutrality. A large lepton 
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Figure 4. lso-5N e g contours in the sin 2 26 s — 8m 2 s plane for L^> — 10 -2 and Sm 2 > (top panel) 
and 8m 2 < (bottom panel), as in Fig. 1. 

asymmetry is responsible for blocking the active-sterile flavor conversions by an in-medium 
suppression of the mixing angle; therefore it has been invoked as a means of significantly 
reducing the sterile abundance [31]. A large lepton number can be generated by e.g. an 
Affleck-Dine mechanism [58] or other models that are able to produce large lepton asymme- 
tries and small baryonic ones [59, 60]. Another interesting possibility is to grow the lepton 
asymmetry from some initial i>) ~ 0(1(T 10 ) using active-sterile oscillations [45, 61, 62]. 
Solving the QKE's in IH and with an initially small but non-zero lepton number, our prelim- 
inary results point toward a final lepton number varying between 10 -5 and 10 -2 depending 



- 11 - 



on the mixing parameters. For illustrative purposes, we choose to adopt = 10~ 2 . 

Figure 4 shows the 5N e s contour plot for = 1CP 2 and 5m 2 s > (top panel) and 
5m? s < (bottom panel). The region with full thermalisation is now much smaller than in 
Fig. 1. 

As we discuss in detail in Appendix A, a large value of confines the resonances to 
very small or large values of x, far away from the maximum of the active neutrino momentum 
distribution (see also [63]). Only at relatively low temperature does the resonance begin to 
move through the momentum distribution. What happens next is qualitatively very different 
for normal and inverted hierarchy. For NH the lepton asymmetry decreases as the resonance 
moves. This causes a run-away effect because as decreases the resonance moves faster, 
causing a faster decrease in L^ a \ When becomes less than approximately 10~ 5 (see 
Eq. (A. 9)), the resonance disappears and the remaining evolution after this point is equivalent 
to the = NH case. For sufficiently large 5m 2 s and sin 2 29 s the non-resonant production 
after the resonance disappears can be significant. However, the required mass difference and 
mixing to obtain the same degree of thermalisation are much larger than in the = 
case. 

The rapid depletion of in NH causes the numerical solution to continue after this 
point with a very small time step. No further resonant production will occur after this point, 
but as we discuss above, some non-resonant thermalisation has yet to happen at this stage. 
To circumvent this problem, we stop the code when becomes very close to zero and 
restart it again with L( fl ) = using the static approximation discussed in Appendix B. 

For IH the lepton asymmetry increases when the resonance moves, causing it to move 
slower and effectively blocking population of the sterile state until very late. For the range 
of mixing parameters studied here production of sterile neutrinos is effectively blocked until 
after the active species decouples, leading to a very small 5N e g. In Appendix A we give 
equations for the position of the resonances for finite along with useful approximations 
valid in different limits. 

For v s mixing parameters as in [56], 5N e s ~ in IH and 5N e s = 0.05 in NH. Constraints 
from BBN, CMB, and LSS have usually assumed a fully thermalised sterile state, but as also 
mentioned in [27] a finite lepton asymmetry can effectively block thermalisation and make 
this assumption invalid. In that case an eV sterile neutrino will not be in conflict with the 
cosmological neutrino mass bound, but of course the extra energy density preferred by CMB 
and LSS will then not be associable with the light sterile neutrino. 

We finally note that since we have solved the quantum kinetic equations using the 1 
sterile + 1 active approximation, only one lepton asymmetry is relevant in our equations 
(either e or fi). However, in the real 3+1 scenario there will be 3 separate flavour asym- 
metries and active-active oscillations will lead to some degree of equilibration between these 
asymmetries. While there will be some quantitative differences between our 1+1 treatment 
and the full 3+1 scenario we do expect the same qualitative behaviour, i.e. a blocking of 
thermalisation due to confinement of the active-sterile resonances. 

4 Conclusions 

Recent cosmological data seem to favor an excess of radiation beyond three neutrino families 
and photons, and light sterile neutrinos are possible candidates. The upcoming measurement 
of 5N e ff by Planck will confirm or rule out the existence of such extra radiation with high 
precision [64, 65]. 
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Light sterile neutrinos could thermalise prior to neutrino decoupling, contributing to 
the relativistic energy density in the early universe. Present data coming from CMB+LSS, 
and BBN allow the existence of one sub-eV mass sterile family but do not prefer extra fully 
thermalised sterile neutrinos in the eV-mass range since they violate the hot dark matter 
limit on the neutrino mass. However, the assumption of full thermalisation is not necessarily 
justified. In this paper, we have studied the evolution of active and sterile neutrinos in the 
early universe in order to calculate the effective number of thermalized species after T ~ 
1 MeV when active neutrinos have decoupled and slightly before BBN commences. We 
have studied the amount of thermalisation for initial null and large (L^ = 10 -2 ) lepton 
asymmetry, for a range of mass-mixing parameters and for both normal and inverted mass 
hierarchies. 

Assuming null initial lepton asymmetry, we find that the assumption of full thermali- 
sation is justified for eV-mass sterile neutrinos with relatively large mixing (as suggested by 
short-baseline oscillation data). This inevitably leads to tension between CMB+LSS data 
which prefers very light sterile neutrinos and Solar, reactor and short-baseline data which 
prefers a mass around 1 eV or higher. 

On the other hand, for large initial lepton asymmetries light sterile neutrinos are not 
(or only partially) thermalised for almost all the scanned parameter space. This provides a 
loophole for eV sterile neutrinos to be compatible with CMB+LSS constraints. For lepton 
asymmetries around 10 -2 almost no thermalisation occurs for the parameters preferred by 
Solar, reactor and short-baseline data, and the sterile neutrinos would contribute very little 
to the current dark matter density. 

One remaining open question neglected in this work is related to the impact of sterile 
neutrinos on BBN. The v e and v e flux distributions are affected by active-sterile conversions 
and they enter the weak rates regulating the neutron-proton equilibrium (see [66] for a review 
on the topic). Therefore the 4 He abundance is sensitive to the presence of sterile families. 
In particular 5N c r > and a less populated v e spectrum are both responsible for increasing 
the freeze-out temperature of the ratio n/p and therefore for a larger 4 He abundance. 

For small and the mixing parameters discussed here the active-sterile oscillations 
occur well before BBN commences, while for large the active-sterile oscillations are no 
longer decoupled from the active ones and can occur close to the BBN temperature. We refer 
the reader to [67] for a discussion of BBN constraints on the sterile sector, but also stress that 
for large values of any quantitative exclusion limits in mixing parameter space would 
require solving the full QKEs including all three active species. This is clearly beyond the 
scope of the present paper, but remains an interesting and important calculation. 

Note added: After the initial version of this paper was finalised, a semi-analytic estimate 
of the BBN effect in the 3+1 scenario using the quantum rate equations has appeared [68]. 
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A Location of the resonances 



Imposing the resonance condition for neutrinos (V z = 0) and for antineutrinos (V z = 0), one 
finds the locations of the resonances [55]. In order to make explicit the x-dependence, we 

define 

Vo 



Introducing 



Vq = — and V\ = V\x 

x 



signfL^] for particles 

— sign[L( a )] for anti-particles 



the resonance conditions (V z = and V z = 0) can be written 

Vix 2 +£\V L \x + V = . 
We define m = sign[5m 2 ] and write the solution in the following way 

x res = x \A£ ± y/X* 



m 



where we have defined 



1 £m 



X 



A 
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(A.3) 
(A.4) 

(A.5) 
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(A.7) 



Note that xo is always real and positive. In order to have a physical solution, F £m has 
to be real and positive. This condition is satisfied for F± t _i(A) for any A and + i(A) for 
A > 1. Thus, we always have two physical solutions when m = — 1, one for particles and one 
for anti-particles. On the other hand, when m = +1 and A > 1, we have two resonances: 
when I > 0, they occur for particles and when £ < 0, they occur for anti-particles being in 
both cases responsible for destroying the lepton number. These equations reproduce the ones 
reported in [55] when m = — 1. 

We can expand the solutions for small and large 
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with A and xq assuming the following expressions 
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For L( fl ) = 10" 2 we have A S> 1 and therefore the lowest resonance will be 
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Note that 

^resjiow is independent on the sign of the mass hierarchy and the total neutrino 
density. Moreover, from the previous equation, we can extract the temperature at which the 
lowest resonance starts sweeping the bulk of the Dermi-Dirac distribution. This provides a 
good estimate of when resonant thermalisation sets in. For example, in the limit of large 
lepton number, assuming x r0Sj i ow ~ 0.1, we find T reSj i ow ~ 3 MeV for (5m 2 , sin 2 29 s ) = 



(1 eV ,10 ). On the other hand, the higher resonance has no effect at all. In fact 



c res,high — Xq X 2A 



180C (3) |l(°)|m. 



77r2T MeV ("^ + n v) 9 



2.6 x 10 



10 



IL^I 
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(n v + n P )gT 2 LeV ' 
Therefore, for Z>) = 0.01, 

^resjhigh w ifi pass through the peak of the Fermi-Dirac distribution 
at T ~ 1 GeV. At that temperature, the damping term is so strong that no oscillations occur 
and thermalisation is inhibited. 



B Adiabatic approximation 

The so-called "adiabatic" approximation was first introduced in [69] and, under certain con- 
ditions, it allows one to derive an approximate analytic solution of the QKE's. In this section, 
we closely follow the derivation from first principles of [49]. Such derivation assumes that 
the rate of repopulation (Pq) vanishes, however the more careful analysis of [70], including 
a non-zero repopulation rate, turns out to give the same final formula for P y . Therefore we 
choose to adopt the simpler derivation. 

Assuming Pq = 0, Eq. (2.4) can be written as a homogeneous matrix equation: 



d_ 
dt 



~Px 








Pz 





-D -V z 
V z ~D -V x 
V~ 



or using a vectorial notation 



dP 



/CP . 



(B.l) 



(B.2) 



The matrix K, can be diagonalised by a time-dependent matrix U, such that lAKlA" 1 = T>. 
The matrix IA defines an instantaneous diagonal basis through Q = UP and, in principle, 
the evolution equation for Q is non-trivial: 



at 



(B.3) 
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However, if we assume that Eq. (B.3) is dominated by the first term, the differential equation 
can be easily solved. This is the so-called "adiabatic" approximation and its applicability 
has been analysed thoroughly in [49]. Quoting [49], it is applicable when 



< 1 



VD 2 + V 2 

T < 3MeV 

cLL^ 



dTMeV 



«5x KT 11 ^ 



MeV 



(B.4a) 
(B.4b) 
(B.4c) 



Equation (B.4a) is not easily stated as just a limit on temperature. If we are not close 
to the resonance and V z is dominated by Vq, we find tan2# s <C 1 which is true for our 
parameter space. If we are close to the resonance, the criterion depends on L^ a ^ through x res 
(see Appendix A for a discussion of the position of the resonances). Using Eqs. (A. 4, A. 8), 
we find 
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(B.5) 
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For large L^ a \ we almost always break the approximation at the lowest resonance. But, 
since the resonance occurs at a very low momentum, it would have no effect on the physics 
anyway. In principle, it could still affect numerics but we did not encounter problems on this 
particular front. For small L", we are safe for most of the parameter space and, as for large 
L( a ), if the resonance is not sitting in a populated part of the Fermi-Dirac distribution, there 
should not be any impact on the physics from breaking this approximation slightly. This 
also applies to the third condition: If the resonance is not in the middle of a populated part 
of the distribution, we do not have a fast evolution of and the approximation is valid. 

Equation (B.3) can be formally solved by 



Qi(t) = exp Ut')dt'^ Qi{t ), 



(B.7) 



where the fcj's are the eigenvalues of K. Expanding those to lowest order in V x , we have 

V 2 D 



h = -D + iV z , 



-D - iV z 



D 2 + V 2 ' 

Assuming that D is large and V x satisfies Eq. (B.4a), we find 

Qi(t) = Q 2 (t) = 0, Q 3 (t) = Q 3 (to). 
The adiabatic approximation allows us to relate P x , P y and P z through 
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where s 3 (t) is the third column in U l {t) which is also the normalised eigenvector corre- 
sponding to A/j. We have 



s 3 (t) = N 



1 

-(D + k 3 )/V z 
-V X (D + h)/(V z k 3 ) 



(B.ll) 



with N a normalisation constant. We can now relate P x and P y to P z to lowest order in V x : 

V X V Z 



Px(t) 



D 2 + V z 2 ~ 

Substituting V z by V z gives the corresponding relations for anti-particles. 



(B.12) 
(B.13) 
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